library(BiocManager)
source("rnaseq_utils.R")
library(DESeq2)
library(BiocParallel)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(tidyr)
library(edgeR)
library(sys)
# remove the outlier
rnaCounts <- rnaCounts %>% select(-`4GEXDARK4`)
sampleAnnotation <- sampleAnnotation[colnames(rnaCounts),]
riboCounts <- riboCounts %>% select(-`4GEXDARK4`)
sampleAnnotation2 <- sampleAnnotation2[colnames(riboCounts),]
rnaCounts = data.frame(rnaCounts)
rownames(sampleAnnotation) = colnames(rnaCounts)
sampleAnnotation$group <- ifelse(substr(sampleAnnotation$group, 1, 1) %in% c("1", "4"),
paste0("X", sampleAnnotation$group),
sampleAnnotation$group)
contrast between day and dark; wild-type
start_time <- Sys.time()
group = factor(sampleAnnotation$group)
y = DGEList(counts = rnaCounts,
group = group,
genes = rownames(rnaCounts))
#y$samples
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes=FALSE]
AveLogCPM <- aveLogCPM(y)
#hist(AveLogCPM)
y <- calcNormFactors(y)
#y$samples
y <- estimateDisp(y, design, robust=TRUE)
#y$samples
#plotBCV(y)
fit <- glmQLFit(y, design, robust=TRUE)
#plotQLDisp(fit)
#summary(fit$df.prior)
day_vs_dark <- makeContrasts(COLENDDAY - COLEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLENDDAY -1*COLEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT5G20630 AT5G20630 5.535501 5.535991 9.724466 7.830552e-15 1.302143e-10
## AT1G69570 AT1G69570 -4.169657 -4.176277 4.964722 5.986341e-14 4.120473e-10
## AT3G09600 AT3G09600 -6.186618 -6.203965 5.354701 8.471142e-14 4.120473e-10
## AT5G05580 AT5G05580 3.479189 3.482670 5.151173 9.911536e-14 4.120473e-10
## AT1G44446 AT1G44446 -3.324531 -3.325194 7.499148 1.560527e-13 4.843196e-10
## AT5G45820 AT5G45820 -6.180542 -6.196619 5.706870 1.747500e-13 4.843196e-10
## AT1G07180 AT1G07180 -4.592792 -4.600768 5.580593 2.696248e-13 6.405131e-10
## AT3G27690 AT3G27690 -4.871827 -4.872516 8.909493 3.335744e-13 6.933760e-10
## AT5G64840 AT5G64840 -3.751972 -3.752280 8.916823 9.196219e-13 1.649972e-09
## AT3G59060 AT3G59060 -3.863711 -3.864380 7.955721 1.012013e-12 1.649972e-09
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLENDDAY -1*COLEXDARK
## Down 1959
## NotSig 12389
## Up 2281
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 4998
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 7.519514 secs
start_time <- Sys.time()
DESeqDataSet = DESeqDataSetFromMatrix(
countData = rnaCounts,
colData = sampleAnnotation,
design = ~ group
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "COLEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 9513
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 951.3
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 4998
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 18.03931 secs
contrast between day and dark; 14G mutant
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X14BENDDAY - X14BEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X14BENDDAY -1*X14BEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT5G45820 AT5G45820 -9.429530 -9.603817 5.706870 5.558380e-14 6.183772e-10
## AT3G27690 AT3G27690 -5.451716 -5.452763 8.909493 8.843228e-14 6.183772e-10
## AT3G47500 AT3G47500 -6.387256 -6.397688 6.374169 1.156934e-13 6.183772e-10
## AT2G05070 AT2G05070 -6.702382 -6.703587 9.838101 1.534871e-13 6.183772e-10
## AT1G69570 AT1G69570 -3.999885 -4.006312 4.964722 1.859334e-13 6.183772e-10
## AT5G57660 AT5G57660 -5.172456 -5.173075 9.479774 3.123604e-13 8.657069e-10
## AT1G18330 AT1G18330 -4.796010 -4.803761 5.709382 4.024328e-13 9.560079e-10
## AT1G44446 AT1G44446 -3.090320 -3.090875 7.499148 5.304359e-13 1.043978e-09
## AT3G59060 AT3G59060 -4.096335 -4.097397 7.955721 5.650253e-13 1.043978e-09
## AT1G07180 AT1G07180 -3.544550 -3.546441 5.580593 2.491965e-12 4.143888e-09
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X14BENDDAY -1*X14BEXDARK
## Down 1880
## NotSig 12443
## Up 2306
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 4954
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.256883 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X14BEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 9667
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 966.7
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 4954
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.470524 secs
contrast between day and dark; 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X4GENDDAY - X4GEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X4GENDDAY -1*X4GEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT1G44446 AT1G44446 -3.953663 -3.954660 7.499148 1.525969e-15 2.537534e-11
## AT5G45820 AT5G45820 -6.767170 -6.793169 5.706870 1.180308e-14 5.976616e-11
## AT1G69570 AT1G69570 -4.069463 -4.075807 4.964722 1.453358e-14 5.976616e-11
## AT3G09600 AT3G09600 -5.679975 -5.691259 5.354701 1.502934e-14 5.976616e-11
## AT1G18330 AT1G18330 -4.701840 -4.706428 5.709382 1.797046e-14 5.976616e-11
## AT3G27690 AT3G27690 -5.016305 -5.017302 8.909493 2.333974e-14 6.468610e-11
## AT3G59060 AT3G59060 -4.007284 -4.007909 7.955721 7.243031e-14 1.720634e-10
## AT3G47500 AT3G47500 -5.279691 -5.286332 6.374169 8.482539e-14 1.763202e-10
## AT5G57660 AT5G57660 -4.701959 -4.702394 9.479774 1.015748e-13 1.876763e-10
## AT2G05070 AT2G05070 -5.289302 -5.289925 9.838101 2.110553e-13 3.266227e-10
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X4GENDDAY -1*X4GEXDARK
## Down 1720
## NotSig 12872
## Up 2037
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 4476
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.352772 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X4GENDDAY", "X4GEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 9236
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 923.6
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 4476
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.394359 secs
contrast between genotypes, day: COL and 14B
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLENDDAY - X14BENDDAY, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLENDDAY -1*X14BENDDAY
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT5G57870 AT5G57870 2.749030 2.749739 7.240689 3.247423e-10 4.619118e-06
## AT5G13210 AT5G13210 -3.824220 -3.838854 3.198047 5.555497e-10 4.619118e-06
## AT5G44610 AT5G44610 2.040782 2.042746 5.585858 1.396575e-08 7.741215e-05
## AT2G18193 AT2G18193 -3.343663 -3.348822 3.958703 5.256919e-08 2.185433e-04
## AT5G60660 AT5G60660 2.159697 2.162756 4.918551 7.241155e-08 2.408263e-04
## AT4G34131 AT4G34131 -1.938129 -1.938531 6.332505 1.258730e-07 3.488570e-04
## AT3G01345 AT3G01345 -3.602524 -3.631246 5.944690 1.924598e-07 4.572020e-04
## AT4G17340 AT4G17340 1.757094 1.757392 7.348337 2.432282e-07 5.055802e-04
## AT4G25210 AT4G25210 -1.541807 -1.542088 6.480979 2.947914e-07 5.221489e-04
## AT3G03920 AT3G03920 -1.887654 -1.887918 6.613104 3.139990e-07 5.221489e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLENDDAY -1*X14BENDDAY
## Down 284
## NotSig 16218
## Up 127
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 623
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.322701 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X14BENDDAY"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 4011
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 401.1
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 623
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.027226 secs
contrast between genotypes, night: COL and 14B
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLEXDARK - X14BEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLEXDARK -1*X14BEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT5G13210 AT5G13210 -3.699243 -3.711515 3.198047 6.465266e-11 1.075109e-06
## AT5G57870 AT5G57870 2.260281 2.260651 7.240689 4.047196e-10 3.365041e-06
## AT5G60660 AT5G60660 2.054796 2.055994 4.918551 2.406564e-09 1.333958e-05
## AT5G14250 AT5G14250 -1.632070 -1.632924 4.981151 3.681860e-09 1.530641e-05
## AT5G66170 AT5G66170 2.316835 2.319725 3.877768 4.723212e-09 1.570846e-05
## AT4G12545 AT4G12545 4.856443 4.889670 3.900695 6.141792e-09 1.702198e-05
## AT4G22490 AT4G22490 2.334567 2.336895 5.016727 1.271862e-08 2.880894e-05
## AT3G21720 AT3G21720 -3.113993 -3.117116 4.479165 1.385961e-08 2.880894e-05
## AT3G11710 AT3G11710 -1.509525 -1.509708 7.236135 2.197153e-08 3.541515e-05
## AT5G43290 AT5G43290 -3.499606 -3.512768 2.504409 2.236003e-08 3.541515e-05
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLEXDARK -1*X14BEXDARK
## Down 289
## NotSig 15992
## Up 348
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 929
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.830286 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X14BEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 5578
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 557.8
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 929
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.36721 secs
contrast between genotypes, day: COL and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLENDDAY - X4GENDDAY, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLENDDAY -1*X4GENDDAY
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT3G01345 AT3G01345 -8.047863 -8.079019 5.944690 3.056858e-14 5.083250e-10
## AT4G15110 AT4G15110 2.878915 2.880247 5.730227 1.877357e-12 1.560928e-08
## AT3G60240 AT3G60240 1.843232 1.843563 7.002984 7.683085e-09 4.258734e-05
## AT1G71030 AT1G71030 2.093515 2.094985 7.033848 2.007862e-07 8.347183e-04
## AT1G80840 AT1G80840 2.046024 2.047474 4.407891 7.730338e-07 2.570956e-03
## AT4G23810 AT4G23810 2.186970 2.190009 3.329671 2.564529e-06 7.107593e-03
## AT1G68840 AT1G68840 1.340636 1.340815 7.553711 6.326135e-06 1.493487e-02
## AT5G03230 AT5G03230 1.488172 1.489980 6.029975 7.184977e-06 1.493487e-02
## AT4G38400 AT4G38400 -1.721493 -1.725036 3.439118 9.668246e-06 1.786370e-02
## AT5G61600 AT5G61600 1.581549 1.582239 5.208417 1.652684e-05 2.748248e-02
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLENDDAY -1*X4GENDDAY
## Down 3
## NotSig 16615
## Up 11
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 15
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.248233 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X4GENDDAY"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 352
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 35.2
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 15
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.037442 secs
contrast between genotypes, night: COL and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLEXDARK - X4GEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLEXDARK -1*X4GEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT3G01345 AT3G01345 -8.971074 -9.050976 5.944690 6.666672e-15 1.108601e-10
## AT3G60240 AT3G60240 1.632912 1.633235 7.002984 1.923606e-07 1.599382e-03
## AT3G10340 AT3G10340 1.638733 1.639960 4.514862 3.316951e-05 1.838586e-01
## AT5G23010 AT5G23010 1.423243 1.423636 7.184509 1.770548e-04 6.469499e-01
## AT3G21720 AT3G21720 -1.718978 -1.721435 4.479165 1.945246e-04 6.469499e-01
## AT5G13210 AT5G13210 -1.626017 -1.635016 3.198047 3.212706e-04 8.904015e-01
## AT2G40330 AT2G40330 -1.978164 -1.983003 3.000185 4.745570e-04 9.779622e-01
## AT2G19970 AT2G19970 -1.405319 -1.406831 3.186981 5.036444e-04 9.779622e-01
## AT5G24240 AT5G24240 -2.511893 -2.532129 1.051515 5.902314e-04 9.779622e-01
## AT2G45360 AT2G45360 -1.608946 -1.615313 1.778893 6.644364e-04 9.779622e-01
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLEXDARK -1*X4GEXDARK
## Down 1
## NotSig 16627
## Up 1
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.288385 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X4GEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 355
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 35.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.407555 secs
contrast between genotypes, day: 14B and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X14BENDDAY - X4GENDDAY, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X14BENDDAY -1*X4GENDDAY
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT4G15110 AT4G15110 2.774434 2.775750 5.730227 4.013744e-12 3.956349e-08
## AT3G01345 AT3G01345 -4.445340 -4.447772 5.944690 4.758373e-12 3.956349e-08
## AT5G57870 AT5G57870 -2.774680 -2.775392 7.240689 1.331251e-10 7.379125e-07
## AT3G60240 AT3G60240 1.830524 1.830854 7.002984 8.904515e-09 3.701829e-05
## AT5G13210 AT5G13210 2.497301 2.502468 3.198047 2.013228e-08 5.781913e-05
## AT3G03920 AT3G03920 2.005648 2.005943 6.613104 2.217572e-08 5.781913e-05
## AT4G12545 AT4G12545 -4.246320 -4.261516 3.900695 2.433904e-08 5.781913e-05
## AT4G34131 AT4G34131 1.860616 1.860990 6.332505 4.188871e-08 8.707092e-05
## AT3G07230 AT3G07230 2.189971 2.190546 5.937643 6.721483e-08 1.098499e-04
## AT2G18193 AT2G18193 2.796600 2.799935 3.958703 7.443747e-08 1.098499e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X14BENDDAY -1*X4GENDDAY
## Down 221
## NotSig 15910
## Up 498
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 1013
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.812748 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X4GENDDAY"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 5629
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 562.9
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 1013
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.369587 secs
contrast between genotypes, night: 14B and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X14BEXDARK - X4GEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X14BEXDARK -1*X4GEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT3G01345 AT3G01345 -4.947873 -4.952555 5.944690 3.957324e-13 6.580635e-09
## AT4G22490 AT4G22490 -2.904688 -2.907204 5.016727 5.414356e-10 4.501766e-06
## AT5G57870 AT5G57870 -2.164870 -2.165234 7.240689 2.806398e-09 1.555587e-05
## AT5G13180 AT5G13180 -2.014076 -2.014778 5.870159 1.080500e-08 4.491907e-05
## AT5G14250 AT5G14250 1.697728 1.698642 4.981151 1.387473e-08 4.614457e-05
## AT3G60240 AT3G60240 1.862560 1.862905 7.002984 2.101704e-08 5.437957e-05
## AT4G12545 AT4G12545 -4.713775 -4.746878 3.900695 2.289115e-08 5.437957e-05
## AT2G31280 AT2G31280 -1.620274 -1.621100 4.937200 4.808899e-08 9.912306e-05
## AT2G45670 AT2G45670 -1.618117 -1.618455 6.162460 5.364770e-08 9.912306e-05
## AT3G07230 AT3G07230 2.405400 2.406257 5.937643 6.682140e-08 1.111173e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X14BEXDARK -1*X4GEXDARK
## Down 289
## NotSig 16062
## Up 278
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 861
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.285785 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BEXDARK", "X4GEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 5456
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 545.6
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 861
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.365066 secs
Ribo-Seq
riboCounts = data.frame(riboCounts)
rownames(sampleAnnotation2) = colnames(riboCounts)
sampleAnnotation2$group <- ifelse(substr(sampleAnnotation2$group, 1, 1) %in% c("1", "4"),
paste0("X", sampleAnnotation2$group),
sampleAnnotation2$group)
group = factor(sampleAnnotation2$group)
y = DGEList(counts = riboCounts,
group = group,
genes = rownames(riboCounts))
y$samples
## group lib.size norm.factors
## X14BENDDAY2 X14BENDDAY 994611 1
## X14BENDDAY4 X14BENDDAY 1128982 1
## X14BEXDARK2 X14BEXDARK 407864 1
## X14BEXDARK3 X14BEXDARK 564737 1
## X14BEXDARK4 X14BEXDARK 332326 1
## X4GENDDAY2 X4GENDDAY 890314 1
## X4GENDDAY3 X4GENDDAY 925586 1
## X4GENDDAY4 X4GENDDAY 1258788 1
## X4GEXDARK2 X4GEXDARK 1151971 1
## X4GEXDARK3 X4GEXDARK 521701 1
## COLENDDAY3 COLENDDAY 234275 1
## COLENDDAY5 COLENDDAY 123746 1
## COLEXDARK2 COLEXDARK 556287 1
## COLEXDARK3 COLEXDARK 2133574 1
## COLEXDARK4 COLEXDARK 350392 1
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes=FALSE]
AveLogCPM <- aveLogCPM(y)
hist(AveLogCPM)
y <- calcNormFactors(y)
y$samples
## group lib.size norm.factors
## X14BENDDAY2 X14BENDDAY 966881 1.0673652
## X14BENDDAY4 X14BENDDAY 1097093 1.0193377
## X14BEXDARK2 X14BEXDARK 398728 0.9595619
## X14BEXDARK3 X14BEXDARK 551606 0.9632180
## X14BEXDARK4 X14BEXDARK 325221 0.9245671
## X4GENDDAY2 X4GENDDAY 861062 1.1586712
## X4GENDDAY3 X4GENDDAY 900301 0.9996082
## X4GENDDAY4 X4GENDDAY 1220563 1.0979268
## X4GEXDARK2 X4GEXDARK 1121272 0.9923085
## X4GEXDARK3 X4GEXDARK 509437 0.8782312
## COLENDDAY3 COLENDDAY 228596 0.9600032
## COLENDDAY5 COLENDDAY 120166 1.1027311
## COLEXDARK2 COLEXDARK 542330 0.9738623
## COLEXDARK3 COLEXDARK 2073598 1.0485569
## COLEXDARK4 COLEXDARK 341882 0.8978050
for (i in 1:ncol(y)){
plotMD(y, column=i)
abline(h=0, col="red", lty=2, lwd=2)
}
contrast between day and dark; wild-type
start_time <- Sys.time()
group = factor(sampleAnnotation2$group)
y = DGEList(counts = riboCounts,
group = group,
genes = rownames(riboCounts))
#y$samples
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes=FALSE]
AveLogCPM <- aveLogCPM(y)
#hist(AveLogCPM)
y <- calcNormFactors(y)
#y$samples
y <- estimateDisp(y, design, robust=TRUE)
#y$samples
#plotBCV(y)
fit <- glmQLFit(y, design, robust=TRUE)
#plotQLDisp(fit)
#summary(fit$df.prior)
day_vs_dark <- makeContrasts(COLENDDAY - COLEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLENDDAY -1*COLEXDARK
## genes logFC unshrunk.logFC logCPM PValue
## AT2G21660 AT2G21660 5.208113 5.216850e+00 8.733563 2.094121e-13
## AT3G63160 AT3G63160 4.949438 4.965324e+00 6.930860 2.415943e-13
## AT3G27690 AT3G27690 -5.388810 -5.396124e+00 9.575708 1.573315e-12
## AT5G20630 AT5G20630 4.556578 4.557859e+00 10.254881 3.134203e-12
## AT3G46970 AT3G46970 5.008629 5.031453e+00 7.569230 3.773531e-12
## AT4G27310 AT4G27310 -11.000459 -1.442695e+08 7.738613 1.020064e-11
## AT5G23240 AT5G23240 6.777450 6.993705e+00 5.901661 1.901113e-11
## AT2G28000 AT2G28000 3.552737 3.555919e+00 8.735217 6.338472e-11
## AT3G45300 AT3G45300 -4.940535 -4.945190e+00 9.635867 6.619262e-11
## AT4G36850 AT4G36850 -5.037495 -5.042810e+00 9.567315 6.718235e-11
## FDR
## AT2G21660 1.336016e-09
## AT3G63160 1.336016e-09
## AT3G27690 5.800287e-09
## AT5G20630 8.347051e-09
## AT3G46970 8.347051e-09
## AT4G27310 1.880318e-08
## AT5G23240 3.003758e-08
## AT2G28000 7.430367e-08
## AT3G45300 7.430367e-08
## AT4G36850 7.430367e-08
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLENDDAY -1*COLEXDARK
## Down 903
## NotSig 9306
## Up 851
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2356
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 4.851323 secs
start_time <- Sys.time()
DESeqDataSet = DESeqDataSetFromMatrix(
countData = riboCounts,
colData = sampleAnnotation2,
design = ~ group
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "COLEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 2977
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 297.7
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2356
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 17.4248 secs
contrast between day and dark; 14B mutant
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X14BENDDAY - X14BEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X14BENDDAY -1*X14BEXDARK
## genes logFC unshrunk.logFC logCPM PValue
## AT4G27310 AT4G27310 -7.611322 -7.731697e+00 7.738613 3.531205e-15
## AT3G47500 AT3G47500 -6.982655 -7.079459e+00 7.313586 2.064323e-14
## AT3G12320 AT3G12320 -7.755058 -7.913509e+00 7.591271 3.447189e-14
## AT3G27690 AT3G27690 -5.296898 -5.302185e+00 9.575708 4.741226e-14
## AT5G06980 AT5G06980 -10.868484 -1.442695e+08 7.519672 4.432385e-13
## AT2G46830 AT2G46830 -8.865830 -9.297847e+00 7.467552 8.195362e-13
## AT3G46970 AT3G46970 5.257768 5.276175e+00 7.569230 1.012280e-12
## AT2G05070 AT2G05070 -5.979803 -5.988081e+00 9.535313 1.865140e-12
## AT5G57660 AT5G57660 -4.473787 -4.477452e+00 9.302224 8.375795e-12
## AT1G06570 AT1G06570 -3.705266 -3.708559e+00 8.484350 4.465806e-11
## FDR
## AT4G27310 3.905513e-11
## AT3G47500 1.141571e-10
## AT3G12320 1.270864e-10
## AT3G27690 1.310949e-10
## AT5G06980 9.804435e-10
## AT2G46830 1.510678e-09
## AT3G46970 1.599402e-09
## AT2G05070 2.578556e-09
## AT5G57660 1.029292e-08
## AT1G06570 4.939182e-08
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X14BENDDAY -1*X14BEXDARK
## Down 870
## NotSig 9258
## Up 932
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2414
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.117298 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X14BEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 3385
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 338.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2414
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.478439 secs
contrast between day and dark; 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X4GENDDAY - X4GEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X4GENDDAY -1*X4GEXDARK
## genes logFC unshrunk.logFC logCPM PValue
## AT4G27310 AT4G27310 -7.701849 -7.806090e+00 7.738613 7.837052e-17
## AT3G12320 AT3G12320 -9.723506 -1.033340e+01 7.591271 3.010892e-16
## AT3G47500 AT3G47500 -7.416789 -7.560539e+00 7.313586 8.625411e-16
## AT5G06980 AT5G06980 -11.222042 -1.442695e+08 7.519672 1.246396e-14
## AT2G46830 AT2G46830 -8.710055 -9.046606e+00 7.467552 2.857466e-14
## AT3G27690 AT3G27690 -4.579005 -4.584166e+00 9.575708 4.131001e-14
## AT1G19530 AT1G19530 -7.344015 -7.425529e+00 7.135410 7.960830e-14
## AT3G16150 AT3G16150 -5.325967 -5.360825e+00 6.352284 4.201473e-13
## AT2G31380 AT2G31380 -5.048744 -5.071614e+00 7.224844 4.328159e-13
## AT1G07040 AT1G07040 -3.911955 -3.920789e+00 7.418560 4.564919e-13
## FDR
## AT4G27310 8.667779e-13
## AT3G12320 1.665024e-12
## AT3G47500 3.179902e-12
## AT5G06980 3.446284e-11
## AT2G46830 6.320715e-11
## AT3G27690 7.614812e-11
## AT1G19530 1.257811e-10
## AT3G16150 5.048800e-10
## AT2G31380 5.048800e-10
## AT1G07040 5.048800e-10
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X4GENDDAY -1*X4GEXDARK
## Down 1235
## NotSig 8588
## Up 1237
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 3210
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.635423 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X4GENDDAY", "X4GEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 4242
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 424.2
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 3210
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.375332 secs
contrast between genotypes, day: COL and 14B
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLENDDAY - X14BENDDAY, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLENDDAY -1*X14BENDDAY
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT5G52310 AT5G52310 2.298437 2.302264 6.242966 3.167691e-06 0.03503466
## AT2G24050 AT2G24050 3.958974 4.001360 5.655656 6.942280e-06 0.03839081
## AT5G11330 AT5G11330 3.469210 3.506837 5.021713 1.431804e-05 0.05278583
## AT5G57870 AT5G57870 2.035819 2.039778 6.949040 3.201044e-05 0.08365051
## AT1G29395 AT1G29395 3.998408 4.026022 5.913189 3.781668e-05 0.08365051
## AT4G36390 AT4G36390 -2.341185 -2.346701 6.177771 7.153424e-05 0.13186144
## AT2G30570 AT2G30570 1.862147 1.862420 10.538731 9.603520e-05 0.13473429
## AT1G51400 AT1G51400 1.689335 1.690210 8.926355 1.112608e-04 0.13473429
## AT1G54410 AT1G54410 1.722246 1.722848 9.235126 1.189421e-04 0.13473429
## AT2G15970 AT2G15970 1.837798 1.838640 8.355813 1.218212e-04 0.13473429
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLENDDAY -1*X14BENDDAY
## Down 0
## NotSig 11058
## Up 2
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 5
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.042515 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X14BENDDAY"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 295
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 29.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 5
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.609078 secs
contrast between genotypes, night: COL and 14B
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLEXDARK - X14BEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLEXDARK -1*X14BEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT4G14130 AT4G14130 3.469643 3.472552 7.859450 1.945077e-09 2.151255e-05
## AT1G55490 AT1G55490 -2.728129 -2.731292 8.973621 3.370916e-08 1.251508e-04
## AT3G56240 AT3G56240 2.332565 2.334391 8.437367 3.924470e-08 1.251508e-04
## AT5G39610 AT5G39610 2.744453 2.751394 6.325387 4.526248e-08 1.251508e-04
## AT1G12900 AT1G12900 -2.433954 -2.434493 9.985552 1.216752e-07 2.691454e-04
## AT3G53460 AT3G53460 -2.755251 -2.759374 8.854840 2.348539e-07 4.329140e-04
## AT3G21720 AT3G21720 -4.232050 -4.291644 4.526621 2.880863e-07 4.551763e-04
## AT4G15390 AT4G15390 3.144799 3.165702 5.863056 5.068291e-07 7.006913e-04
## AT1G74880 AT1G74880 -3.140100 -3.176221 5.634535 9.189849e-07 9.218939e-04
## AT3G46010 AT3G46010 2.266082 2.273641 6.647010 9.443118e-07 9.218939e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLEXDARK -1*X14BEXDARK
## Down 193
## NotSig 10710
## Up 157
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 589
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.000142 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X14BEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 1786
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 178.6
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 589
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.342794 secs
contrast between genotypes, day: COL and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLENDDAY - X4GENDDAY, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLENDDAY -1*X4GENDDAY
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT3G60240 AT3G60240 2.189418 2.192475e+00 7.754499 3.920502e-07 0.004336075
## AT4G36850 AT4G36850 -2.817004 -2.821714e+00 9.567315 3.956657e-06 0.021880315
## AT3G15356 AT3G15356 -2.593780 -2.600656e+00 7.013375 6.057575e-05 0.223322582
## AT5G20150 AT5G20150 2.810670 2.831369e+00 4.183695 1.308052e-04 0.301615767
## AT1G02640 AT1G02640 -1.951093 -1.954489e+00 8.492358 1.363543e-04 0.301615767
## AT2G17230 AT2G17230 -1.966528 -1.970363e+00 6.563663 2.130801e-04 0.344802627
## AT2G15080 AT2G15080 -7.394377 -1.442695e+08 4.578464 2.182295e-04 0.344802627
## AT3G45970 AT3G45970 -2.298422 -2.304392e+00 6.588248 2.583441e-04 0.357160767
## AT1G22770 AT1G22770 2.720781 2.751204e+00 4.236109 3.218990e-04 0.395578152
## AT1G15430 AT1G15430 -7.127060 -1.442695e+08 4.353627 5.250017e-04 0.548527832
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLENDDAY -1*X4GENDDAY
## Down 1
## NotSig 11058
## Up 1
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.044354 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X4GENDDAY"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 91
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 9.1
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.016833 secs
contrast between genotypes, night: COL and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(COLEXDARK - X4GEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*COLEXDARK -1*X4GEXDARK
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT3G60240 AT3G60240 1.862816 1.864929 7.754499 5.438739e-06 0.06015245
## AT4G20390 AT4G20390 3.323862 3.403996 4.126823 1.316746e-04 0.60463238
## AT3G27160 AT3G27160 1.680439 1.684490 7.616775 2.323116e-04 0.60463238
## AT2G40330 AT2G40330 -2.192904 -2.206413 5.029629 2.573583e-04 0.60463238
## AT3G25730 AT3G25730 -3.094740 -3.195454 3.589990 2.733419e-04 0.60463238
## AT2G43920 AT2G43920 -1.822146 -1.831592 4.344406 4.399068e-04 0.75017639
## AT5G09220 AT5G09220 -2.151239 -2.169488 4.732271 4.747952e-04 0.75017639
## AT5G23010 AT5G23010 1.666882 1.671034 7.404208 1.001170e-03 1.00000000
## AT2G23170 AT2G23170 -1.651932 -1.655450 5.954169 1.103876e-03 1.00000000
## AT4G14130 AT4G14130 1.487309 1.487828 7.859450 1.195258e-03 1.00000000
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*COLEXDARK -1*X4GEXDARK
## Down 0
## NotSig 11060
## Up 0
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 1
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.103256 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X4GEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 25
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 2.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 1
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.344416 secs
contrast between genotypes, day: 14B and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X14BENDDAY - X4GENDDAY, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X14BENDDAY -1*X4GENDDAY
## genes logFC unshrunk.logFC logCPM PValue FDR
## AT3G60240 AT3G60240 2.121398 2.124404 7.754499 1.508249e-07 0.001409122
## AT5G57870 AT5G57870 -2.368533 -2.372729 6.949040 3.711322e-07 0.001409122
## AT2G24050 AT2G24050 -3.950367 -3.992707 5.655656 4.412224e-07 0.001409122
## AT5G17710 AT5G17710 2.269856 2.273933 6.150627 5.096283e-07 0.001409122
## AT2G18193 AT2G18193 3.057497 3.101369 4.288943 1.230352e-06 0.002721539
## AT5G36910 AT5G36910 -2.650604 -2.662379 7.378530 2.130043e-06 0.003926379
## AT4G15110 AT4G15110 5.222394 5.324235 4.877607 4.138304e-06 0.006515614
## AT5G11330 AT5G11330 -3.048942 -3.085444 5.021713 4.712922e-06 0.006515614
## AT2G22360 AT2G22360 1.838789 1.845057 5.169158 2.214719e-05 0.027216437
## AT5G08610 AT5G08610 2.370609 2.381078 4.951388 2.869953e-05 0.030705817
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X14BENDDAY -1*X4GENDDAY
## Down 7
## NotSig 11043
## Up 10
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 40
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.072463 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X4GENDDAY"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 502
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 50.2
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 40
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.390091 secs
contrast between genotypes, night: 14B and 4G
start_time <- Sys.time()
day_vs_dark <- makeContrasts(X14BEXDARK - X4GEXDARK, levels = design)
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient: 1*X14BEXDARK -1*X4GEXDARK
## genes logFC unshrunk.logFC logCPM PValue
## AT3G56240 AT3G56240 -2.770708 -2.772654 8.437367 3.163044e-09
## AT3G12780 AT3G12780 2.651389 2.651975 10.702835 3.073265e-08
## AT3G27160 AT3G27160 2.761626 2.766593 7.616775 1.052553e-07
## AT1G55490 AT1G55490 2.874402 2.877830 8.973621 1.329417e-07
## AT2G43030 AT2G43030 2.204186 2.205813 8.520207 4.697147e-07
## AT5G39610 AT5G39610 -2.537727 -2.544460 6.325387 5.685509e-07
## AT3G51920 AT3G51920 -1.991011 -1.992451 7.715316 5.728854e-07
## AT3G46010 AT3G46010 -2.438612 -2.446364 6.647010 6.317070e-07
## AT3G24430 AT3G24430 3.012059 3.028821 6.325834 7.374642e-07
## AT2G41680 AT2G41680 2.096484 2.099069 7.449179 1.107221e-06
## FDR
## AT3G56240 3.498327e-05
## AT3G12780 1.699516e-04
## AT3G27160 3.675838e-04
## AT1G55490 3.675838e-04
## AT2G43030 8.733350e-04
## AT5G39610 8.733350e-04
## AT3G51920 8.733350e-04
## AT3G46010 8.733350e-04
## AT3G24430 9.062616e-04
## AT2G41680 1.054867e-03
is.de <- decideTestsDGE(tr)
summary(is.de)
## 1*X14BEXDARK -1*X4GEXDARK
## Down 133
## NotSig 10790
## Up 137
plotMD(tr, status=is.de)
output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 459
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.07013 secs
start_time <- Sys.time()
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BEXDARK", "X4GEXDARK"))
# indexes that all have a value
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 1425
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 142.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 459
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.362078 secs
